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Quantum optimal control theory is a powerful tool for engineering quantum sys- 
tems subject to external fields such as the ones created by intense lasers. The 
Qh' formulation relies on a suitable definition for a target functional, that translates the 

Slj ' intended physical objective to a mathematical form. We propose the use of target 

^ . , functionals defined in terms of the one-particle density and its current. A strong 

motivation for this is the possibility of using time-dependent density-functional the- 

cn 

J> ' ory for the description of the system dynamics. We exemplify this idea by defining 

o^ : 

CN ■ an objective functional that on one hand attempts a large overlap with a target 

density and on the other hand minimizes the current. The latter requirement leads 
to optimized states with increased stability, which we prove with a few examples of 
one- and two-dimensional one-electron systems. 
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I. INTRODUCTION 



Quantum control [l| is concerned with the detailed manipulation of systems in the micro- 
scopic world. Its most common application is the design of laser fields capable of triggering a 
given response of some target system. Under this umbrella name of quantum control we may 
include both experimental techniques, the theoretical description of those, and the mathe- 
matical techniques used to predict the controlling external fields. Experimentally, the field 
has advanced very rapidly in the last years due to the convergence of several important paths: 



most notably, the_deyelopment of ultra-fast and ultra-intense laser sources, t 
of pulse shapers 



le appearance 



12|. 



-7|, and the invention of the adaptive feedback technique 
Theoretically, the most general framework used to address quantum control problems 
is quantum optimal control theory (QOCT) 13l-ll6l|. It is formulated as the problem of 
maximizing a functional that quantifies to what extent the given objective is met by a 
certain laser field. The result is typically an algorithm that requires, as main steps, the 
propagation of the time dependent Schrodinger equation that describes the evolution of the 
system (and, usually, of a similar equation that is propagated backwards). The ab initio 
solution of these equations is, unfortunately, prohibitively expensive for most systems due 
to the complexity of their wavefunctions. 

Time-dependent density-functional theory (TDDFT) 



13, 



18| provides in many cases a 



viable solution for this problem. Recently, the merger of QOCT and TDDFT has been 
proposed and numerically demonstrated 19|] . Within TDDFT the electron density, that can 
be obtained from the solution of the so-called time dependent Kohn-Sham (KS) equations, 
substitutes the many-electron wavefunction as the main object to be manipulated. 

The targets of QOCT are usually defined as functionals of many-electron wavefunctions. 
For example, if one wishes to increase the population of a given excited state, the target is 
defined in terms of the projector onto that particular excited state wave function. However, 
these wavefunctions are in principle unavailable within the TDDFT framework that provides 



the electron density only 



201 ]. Therefore, if one is to use QOCT with TDDFT consistently. 



the targets should be defined as functionals of this electron density alone. 

Yet, to construct better target functionals (whether or not TDDFT is used), it may be 
useful to add an extra ingredient to the electron density, namely the electron current density. 
Note that the longitudinal component of this current vector field is linked directly to the 
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density through the continuity equation, and can therefore be obtained with TDDFT (the 
perpendicular component of the KS current density is also conjectured to coincide with the 
many-electron one). 

For example, it is easy to think of a density-functional target constructed to maximize 
charge transfer between two regions of a system. However, once the controlling field is 
switched off, even if the charge has been effectively transferred, the final state may not be 
stable. The inclusion of the current in the definition of the target functional can be used, 
as will be demonstrated below, to stabilize the final state of the propagation: if the current 
(or its divergence) is small, the continuity equation ensures that the time-variation of the 
density is also small. 

In order to demonstrate the usefulness and feasibility of including the electron current 
density in the definition of target functionals, we studied several model one-electron systems 
in one and two dimensions. The point is to show that with this combined target one can 
obtain solutions that (i) are stabler if the current is minimized and that (ii) resemble results 



as if the target was formulated in terms of the wavefunction [2]J]. Finally, we also note 
that we do not put forward the use of current TDDFT 22], but rather the inclusion of the 
current density in the definition of the control targets. The proposed equations have been 
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2^. 



implemented into the octopus code 

The paper is organized as follows: in Section |Tl] we briefly outline the optimization scheme, 
and introduce the new control functionals defined in terms of the electron density and its 
current. In Section IIIII we present numerical examples that serve as a proof of principle. 
The conclusions are given in Section |IVl Atomic units (a.u.) are used throughout. 



II. FORMALISM 

We consider a one-electron system governed at rest by the Hamiltonian Hq that interacts 
with a time-dependent control field e{t) in the time span t G [0,T], such that the total 
Hamiltonian is given by: 

H{t) = Ho- fMe{t), 0<t<T, (1) 

where /i is the dipole operator that couples with the electric field e(t). Without loss of 
generality, we will think of e{t) as the electric field of a laser pulse. 
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In QOCT one tries to find the control field that best realizes an objective that in general 
is mathematically encoded into a functional of the field G[e\. This is defined in the following 
way: the field determines the evolution of the system, e — ^[e], and this in turn determines 
G through an intermediate functional F: 



G[e] 



e ,e 



(2) 



The functional F must be carefully defined in order to ensure the fulfillment of the objective. 
Typically, it is split into two pieces: 



(3) 



where 0[%l)\ is the quantity that we really want to optimize, whereas J-" imposes a penalty on 
undesired features of the controlling fields: e.g. high frequency components, unrealistically 
high intensities, etc. In our case we have used the following expression: 

T 

J^[e] = - j dta{t)e^{t). (4) 



If a = 1, J-" would be minus the fluence, or integrated energy of the laser pulse. Therefore, 
this term penalizes pulses with too high intensities. The function a may then be chosen 
to ensure that the laser pulse is smoothly switched on and off, by penalizing non-zero field 
values near the beginning and the end of the pulse: 



<|erf 
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(5) 



where erf is the error function. 

The search for the maxima of G is substituted by the search for the maxima of F, which 
however cannot be unconstrained, since the evolution of the system must obey Schrodinger's 
equation. The formalism takes care of this by introducing a Lagrange functional: 

T 



-2^He 



{x{t)\^^+iHo-ifie{t)m)) 



(6) 



where the auxiliary state x is introduced as a Lagrange multiplier. The functional whose 
critical points are to be found is: 



(7) 
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The corresponding Euler-Lagrange equations are: 

d 



t^Mt) = [H^- ^ie{t)]ip{t), (8a) 



m = V^o, (8b) 

'dr 



i^Y{t) = [Ho - I2e{t)] x{t) , (8c) 



ait)e{t) = -3m (xitMm)- (8e) 

Equation (l8a|) is merely Schrodinger's equation for the system under the influence of 
the controlling laser e. It must be solved by forward propagation, since it is an initial value 
problem with boundary condition given by Eq. (ISbp . Equation (IHcj) is also a Schrodinger-like 
equation for the Lagrange multiplier state x- However, in this case it must be propagated 
backwards since the boundary condition is given at the final propagation time T. Finally, 
Eq. ( I8e|) couples both t/j and x, and permits to obtain the solution laser field e(t) (the dipole 
matrix element between states x ^"^^ "0 is the principal ingredient). 

These equations are coupled, and their self-consistent solution provides the critical points 
of the functional G. In order to find one solution that corresponds to a maximum, we have 
in this work utilized the monotonic algorithm described in Ref. 
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It consists of a series of 



back- and forward propagations. For other possible algorithms to solve Eqs. ( IHal) - ( IHel) we 



refer the reader, for example, to Refs. 



26 



and 
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Now we focus our attention to the precise form of the functional 0[iIj] (generally speaking, 
it is a functional of the full evolution of the wave function i/j; however in the previous 
equations we have assumed it depends only on the final state ip{T)). As discussed in the 
introduction above, we will assume that it is a functional of the density n and of the current 
density j: 

0[i;] = 0[n[i;{T)],j[i;{T)]], (9) 
where n and j, for the one-electron case, read: 

n[tp]{r) ='^*{r)tp{r) , (10a) 
j[i(j]{r) = 3m [ij*{r)V'ilj (r)] . (10b) 

The particular form for O given in Eq. (|9]) leads to the following expression for the final- value 
condition for the auxiliary wave function x [Eq. ( l8d|l ]: 
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5tlj*{r,T) 



SO 



Sn{r,T) 



Viljir,T) + -ijj{r,T)V 



60 



5jir,T) 



, (11) 



This expression is valid for any target functional defined in terms of the density and its 
current. Let us narrow it down for a particular case: we wish to maximize the overlap of 
the density n with a target density ritg (at the end of the propagation). We may define for 
that purpose a first functional Oi. 

2 

^|n{r)- ^ln,Jr) . (12) 



It ranges in the interval [—2, 0] with its maximum at complete overlap. However, if we only 
use this functional, the resulting density at the end of the propagation will not be stable and 
change rapidly after the controlling field is switched off. In order to stabilize the achieved 
state, we may add a current dependent functional, making use of the fact that the temporal 
behavior of the density is connected with the current j{r,t) by the continuity equation: 

dn{r, t) 



dt 



(13) 



For any eigenstate of Ho the divergence of the current vanishes, and the density remains 
constant. Suppressing the current will assure a stationary density, at least right after the 
controlling pulse is switched off, and perhaps will ensure small variations thereafter. There- 
fore we implemented an objective functional whose maximum is at zero current, namely: 



OM = -Wc jd?r\j{r)\\ 



(14) 



Here, Wc > is weighting the importance of the current suppression, which is useful if one 
combines the objectives in Eqs. (fT2l) and (fT4|) : 



0[^] = 0[n[^{T)imT)]] 

= OMi^iT)]] + o,[mT)]] 



(15) 



In this case Wc should also homogenize the dimensions of the functionals. It only remains 
to rewrite Eq. ffTTj) for this particular case: 



50 



n{r,T) 



2iwJ{r, T) ■ V^(r, T) + iw,i:{r, T) V • j{r, T) . (16) 
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Two remarks are in order before moving on to applications: First, instead of minimizing 
the current, we also applied this scheme to minimizing the divergence of the current, since 
this is the quantity to which the temporal behavior of n is directly related by Eq. (fT3l) . 
However, the more stringent minimization of the total current turned out to be numerically 
preferable. Finally, note that in the one-particle case, the wave function ifj can be written as 
il){r,t) = ^Jn{r, t)e^^^'^'^\ and then the current density is given by j{r^t) = n{r , t)'V S {r , t) . 
Hence controlling the density and the current density amounts to controlling the full wave 
function directly, since both objects contain the same information. 



III. APPLICATIONS 



In the following, we present three illustrative examples of the method described above: 
electron transfer in both a ID and a 2D asymmetric double well, and the Is — transition 
in a 2D model of the hydrogen atom. In all calculations, there are a few parameters and 
settings whose values have an influence in the final outcome, most notably the value of the 
weight Wc and the initial guess for the laser field. Regarding the former, we observed that 
a value of around Wc = 10 led to reasonable convergence and stable results. Regarding the 
initial guess, we found it advisable to do several calculations with different starting points, 
since the search space contains many local minima. 



A. ID asymmetric quantum well 

We consider a ID asymmetric double well, defined by the following potential function: 



which is sometimes used in quantum chemistry to model isomerisation processes 28|. Note 
that even if we are usually referring to electronic processes and electronic densities and 
currents, it need not be the case, and the wave function may be that of a nuclear wave 
packet. The goal is to guide a transition from the ground state (GS) to a state whose 
density is equal to that of a superposition of the GS and the first excited state (ES), within 
a total time of T = 300 a.u. We choose the superposition, 

^tg(x) = ^(^GS(X) + i^Esix)), (18) 
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FIG. 1. (Color online) ID QW. The overlap between the controlled density n and the target density 
ntg is shown in the top panel. For minimized current {wc = 10) the overlap with the target (found 
in the inset) at terminal time T = 300 a.u. is better. Additionally, density n is significantly more 
stable than without current control {wc = 0). The respective laser pulses are reported in the lower 
panel. 



whose corresponding density ntg = | iptg \ ^ is illustrated in the inset of the top panel of Fig. [H 

As the target density ritg does not correspond to that of a stationary state, it is crucial 
not only to maximize the overlap according to Eq. (IT^ . but also to minimize the current 
by setting Wc = 10 in Eq. fll4p . The top panel of Fig. [T] shows for the cases of Wc = 10 and 
Wc = that a very good overlap is found at final time T = 300 a.u. In fact, mapping the 
target functional values on the interval [0, 1], we obtain 0.99998 for the first and 0.99989 for 
the second case. However, minimizing the current guarantees a significantly longer lifetime 
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of the density in its target shape. The two respective laser fields are found in the bottom 
panel. 



B. 2D asymmetric quantum well 

A 2D asymmetric double quantum well can be realized by adding a parabola in the 
?/-direction to the previous ID potential: 

V(x, y) = —x^ + —x^ - -x^ + -y^ . (19) 
^ '^^ 64 256 4 2^ ^ ^ 

The top panel of Fig. [2] displays the densities of the GS and the first ES (indicated as ntg), 

both as ID plots along the x-axis {y = 0) and as density plots in the xy-plane (in the insets). 

The objective of this example is to perform a density transfer from the GS to the ES with 

and without suppressing the current, i.e. with Wc = 20 and Wc = 0, respectively. 

As can be see in the top panel of Fig. |2] both optimizations lead to a very good overlap 

with the target density ritg at T = 300 a.u. However, the intensities of the optimizing laser 

fields were rather different: in the case where the minimization of the current was enforced, 

the electric field was roughly half the size of its counterpart without current control. Also, 

the behavior of the respective densities diverges considerably once the laser is switched off. 

This is illustrated in Fig. [3] where we show the difference between the target ritg and the 

controlled density n at two crucial points along the x-axis. While the top panel does so for 

the node of ritg at x = —1.5 a.u., the bottom panel refers to its maximum at x = 2.4 a.u. The 

requisite of current suppression induces much smaller oscillations of the controlled density 

around a value closer to the target. 



C. 2D hydrogen: transition Is — )■ 2pj^ 



Let us now consider a soft-Coulomb 2D model of the hydrogen atom. In this case the 
objective is a stable density after a transfer from the Is GS to the 2px ES. The respective 
soft-Coulomb potential reads: 

-1 



V{x,y) 



(20) 



^yl + x^ + y^ 

and is shown in the bottom panel of Fig. |H Such soft potentials are used to model Wannier 



excitons 



to describe multi-photon processes 



30| and to understand strong laser fields 
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FIG. 2. (Color online) 2D QW. The top panel shows the optimized densities (with and without 
current suppression) as well as the densities of the GS and of the target ntg along the x-axis 
(y = a.u.). Both optimizations achieve an excellent target density overlap. The insets show the 
respective 2D plots of the GS and ES density. The bottom panel shows the optimal lasers. 



31| . The densities of the initial and target state ntg are shown in the top panel of 



better 
Fig.H 

At final time T = 700 a.u. the overlap of the controlled density n with the target density 
TT-tg is similar between the optimization with current suppression [wc = 40) and without 
{wc = 0) as shown in the top panel of Fig. O It is good for x > a.u. (especially if current 
suppression enforces a solution with a node at a; = a.u.), while it remains unsatisfactory 
for X < a.u. We underline that this asymmetry in attaining the target is a consequence of 
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FIG. 3. (Color online) 2D QW. The difference between the controlled density n and the target 
density ntg as a function of time at the node (x = —1.5 a.u., top panel) and at the maximum 
(x = 2.4 a.u., bottom panel) of the target. The colored background marks the time t > 300 a.u. 
without laser. The insets in both panels are guides to the eye. They show the target density with 
a mark on the node and on the maximum, respectively. 

persisting density oscillations as explained below. Accordingly, the region of good overlap is 
periodically changing from on side of the x-axis to the other. Nonetheless, this asymmetry is 
influenced by the initial guess laser field. In fact, conserving the frequency of the guess field, 
but changing the sign of its amplitude leads to the optimal laser with inverted sign. This 
"negative" optimal laser causes in turn exactly the opposite terminal configuration: with a 
good overlap for x < a.u. and poor overlap otherwise (not shown). Note that in terms 
of the objective functional the final density configuration and its "mirrored" counterpart 
are equally valid. The respective optimal lasers for Wc = 40 and Wc = are plotted in the 
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bottom panel. 

However, a substantial difference between the two optimal solutions can be seen in their 
behavior once the control is switched off. Although we observe an oscillation of density 
between the two lobes in both cases, the minimization of the current leads to a significantly 
smaller rate of fluctuation. Figure [6] illustrates the difference between the controlled and 
target density at one maximum (at x = —2 a.u., top panel) and at the node (at x = a.u., 
bottom panel) of Utg. The difference between the target and the controlled density at the 
node is reduced by a factor of five in the case of minimized current. This means that this 
point of zero density is well represented if the current is held small. Consequently, only 
a small portion of density can pass from one side to the other. Indeed, oscillations of the 
difference at the maximum are halved if the current is suppressed. 

A (soft-) Coulomb potential implies the presence of degeneracies and close-lying states, 
e.g. the target state is degenerate in energy with 2py, and a large number of energetically 
close states is available. Therefore it is impossible to totally exclude a quantum mechanical 
superposition of states (since the target is not a given state, but merely a given density), and 
this is the origin of the persisting oscillations in the optimal density. It is interesting to note 
that even our additional calculations that involved a target in terms of the wavefunction 21 1 
lead to similar optimal solutions with the same problem. Finally, we also stress that this 
kind of potential is fundamentally different from the potentials used in the quantum well 
examples, since it vanishes asymptotically (and does not grow to infinity) at large distances. 
Numerically, this fact may imply the appearance of undesired border effects, and therefore 
we had to use a much larger radius for the simulation box (r = 160 a.u.). 



IV. CONCLUSIONS 

We have proposed the use of target functionals defined in terms of both the one-particle 
density and the current density for QOCT calculations. In particular, we have shown the 
use of functionals that maximize the overlap of the controlled density with a given target 
density, and simultaneously minimizes the current, in order to stabilize the final state at 
the end of the action of the controlling laser pulse. Such an objective fits very well into a 
TDDFT description of the system. A proof of concept was offered with three prototypical 
ID and 2D systems. In these cases, we observed how the suppression of the current reduced 
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FIG. 4. (Color online) 2D H-atom. The top panel illustrates the densities of the Is GS and the 
2px ES along the x-axis (y = a.u.), their 2D plots are shown in the inset. The target state has a 
node at rr = a.u. The bottom panel shows the soft-Coulomb potential V along the x-axis (y = 0). 

quantum oscillations of the final state. 
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